﻿#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <vector>

using namespace std;

/**
 * @brief 辅助函数: 从 CSV 文件加载数据
 *
 * @param path
 * @param delimiter
 * @return vector<vector<double>>
 */
vector<vector<double>> loadCSV(const string &path, char delimiter = ',')
{
    ifstream file(path);
    if (!file)
        throw runtime_error("Cannot open file: " + path);
    string                 line;
    vector<vector<double>> data;
    while (getline(file, line))
    {
        stringstream   ss(line);
        string         cell;
        vector<double> row;
        while (getline(ss, cell, delimiter))
        {
            if (!cell.empty())
                row.push_back(stod(cell));
        }
        if (!row.empty())
            data.push_back(row);
    }
    return data;
}

/**
 * @brief
 *
 * @param X
 * @return vector<double>
 */
vector<double> mean(const vector<vector<double>> &X)
{
    int            n = X.size(), m = X[0].size();
    vector<double> mu(m, 0.0);
    for (auto &row : X)
        for (int j = 0; j < m; ++j)
            mu[j] += row[j];
    for (double &v : mu)
        v /= n;
    return mu;
}

vector<vector<double>> center(const vector<vector<double>> &X, const vector<double> &mu)
{
    int                    n = X.size(), m = X[0].size();
    vector<vector<double>> Y = X;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            Y[i][j] -= mu[j];
    return Y;
}

vector<vector<double>> covariance(const vector<vector<double>> &X)
{
    int                    n = X.size(), m = X[0].size();
    vector<vector<double>> S(m, vector<double>(m, 0.0));
    for (int i = 0; i < n; ++i)
        for (int a = 0; a < m; ++a)
            for (int b = 0; b < m; ++b)
                S[a][b] += X[i][a] * X[i][b];
    double scale = 1.0 / (n - 1);
    for (int a = 0; a < m; ++a)
        for (int b = 0; b < m; ++b)
            S[a][b] *= scale;
    return S;
}

/**
 * @brief 向量点积与范数
 *
 * @param a
 * @param b
 * @return double
 */
double dot(const vector<double> &a, const vector<double> &b)
{
    double s = 0;
    for (size_t i = 0; i < a.size(); ++i)
        s += a[i] * b[i];
    return s;
}

/**
 * @brief 标准化
 *
 * @param a
 * @return double
 */
double norm(const vector<double> &a)
{
    return sqrt(dot(a, a));
}

/**
 * @brief 矩阵向量乘法
 *
 * @param A
 * @param x
 * @return vector<double>
 */
vector<double> matvec(const vector<vector<double>> &A, const vector<double> &x)
{
    int n = A.size(), m = A[0].size();

    vector<double> y(n, 0.0);

    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            y[i] += A[i][j] * x[j];
    return y;
}

/**
 * @brief 幂迭代求特征向量
 *
 * @param A
 * @param max_iter
 * @param eps
 * @return pair<double, vector<double>>
 */
pair<double, vector<double>> powerIteration(const vector<vector<double>> &A, int max_iter = 1000, double eps = 1e-6)
{
    int            n = A.size();
    vector<double> b(n, 1.0);
    double         lambda = 0;
    for (int it = 0; it < max_iter; ++it)
    {
        vector<double> b_next = matvec(A, b);
        double         norm_b = norm(b_next);
        for (double &v : b_next)
            v /= norm_b;
        double lambda_next = dot(b_next, matvec(A, b_next));
        if (fabs(lambda_next - lambda) < eps)
            break;
        b = b_next;
        lambda = lambda_next;
    }
    return {lambda, b};
}

/**
 * @brief 对称矩阵的 deflation（去掉已求特征）
 *
 * @param A
 * @param v
 * @param lambda
 */
void deflate(vector<vector<double>> &A, const vector<double> &v, double lambda)
{
    int n = A.size();
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            A[i][j] -= lambda * v[i] * v[j];
}

/**
 * @brief PCA 主函数
 *
 * @param X
 * @param k
 * @return vector<vector<double>>
 */
vector<vector<double>> pca(const vector<vector<double>> &X, int k)
{
    int                    n = X.size(), m = X[0].size();
    vector<double>         mu = mean(X);
    vector<vector<double>> Xc = center(X, mu);
    vector<vector<double>> S = covariance(Xc);

    vector<vector<double>> components;  // 每列为主成分
    for (int i = 0; i < k; ++i)
    {
        auto [lambda, v] = powerIteration(S);
        components.push_back(v);
        deflate(S, v, lambda);
    }

    // 投影
    vector<vector<double>> Z(n, vector<double>(k, 0.0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < k; ++j)
            for (int t = 0; t < m; ++t)
                Z[i][j] += Xc[i][t] * components[j][t];
    return Z;
}

/**
 * @brief 主函数
 *
 * @return int
 */
int main()
{
    try
    {
        auto X = loadCSV("iris.csv");
        cout << "数据维度: " << X.size() << "x" << X[0].size() << endl;
        cout << "前5个样本:" << endl;
        for (int i = 0; i < 5; ++i)
        {
            for (double v : X[i])
                cout << v << " ";
            cout << endl;
        }

        auto Z = pca(X, 2);
        cout << "降维后前5个样本:" << endl;
        for (int i = 0; i < 5; ++i)
        {
            for (double v : Z[i])
                cout << v << " ";
            cout << endl;
        }
    }
    catch (const exception &e)
    {
        cerr << "Error: " << e.what() << endl;
        return 1;
    }
    return 0;
}
